---
title: "04_Seurat_DA"
author: "Jae"
date: "2023-06-27"
output: html_document
---

# Load libraries
```{r}
library(here)
library(Seurat)
library(tidyverse)
library(biomaRt)
library(dplyr)
library(gridExtra)
library(patchwork)
library(ggplot2)
library(grid)
library(data.table)
library(RColorBrewer)
library(monocle3)
library(Matrix)

"config.r" %>%
    here() %>%
    source()

cr.dir <- file.path(cellranger_multi_out_dir, library_name, "outs", "per_sample_outs")

```

# Load Seurat Objects
```{r}
# load in all seurat object
in.dir <- "~/Desktop/Bioinformatics_CryptoBrain/Stefs_scRNA-seq_data/Updates_Jae"
setwd(in.dir)


seu.file <- "seu_obj_cc_020924.rds"
seu.obj.cc <- readRDS(seu.file) 
```


# Metadata Update (seu.obj.cc)
```{r}
# Update metadata table for **seu.obj.cc**

# split "orig.ident" into column variable "treatment" (i.e. WT, GENE-X) 
# and column variable "time.point" (e.g. 0d, 4d, 14d)
seu.obj.cc@meta.data <- seu.obj.cc@meta.data %>%
  separate(col = orig.ident, into = c("treatment", "time.point"), sep = "_")

# Replace "no" label with "0d" under time.point
seu.obj.cc@meta.data <- seu.obj.cc@meta.data %>%
  mutate(time.point = replace(as.character(time.point), time.point == "no", "0d"))
```


# Construct UMAPs (seu.obj.cc)

seu.obj.cc
```{r}
# n = number of Cells
n = 16000

seu.obj.cc <- FindVariableFeatures(seu.obj.cc, selection.method = "vst")
seu.obj.cc <- ScaleData(seu.obj.cc, selection.method = "vst")
seu.obj.cc <- RunPCA(seu.obj.cc, verbose = FALSE)

seu.obj.cc <- Seurat::FindNeighbors(seu.obj.cc, 
                                    k.param = sqrt(n), 
                                    dims = 1:30
                                    )

seu.obj.cc <- Seurat::FindClusters(seu.obj.cc,
                                   algorithm = 1,
                                   resolution = 0.7,
                                   random.seed = 72423
                                   )

seu.obj.cc <- RunUMAP(seu.obj.cc, reduction = "pca", dims = 1:30)

Seurat::DefaultAssay(seu.obj.cc) <- "RNA"
Idents(seu.obj.cc) <- "seurat_clusters"


seu.obj.cc.plot <- DimPlot(seu.obj.cc, label = TRUE, label.size = 4.5) + 
  theme(legend.position = "none") + 
  ggtitle("CC-Pooled-UMAP") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-15,20), ylim = c(-15,20))

feat.plot <- FeaturePlot(seu.obj.cc, features = "Spp1") + 
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-10,20))

grid.arrange(seu.obj.cc.plot, feat.plot, ncol=2)
```



# Subset pooled integrated seu object (seu.obj.cc) and plot
```{r}
Idents(seu.obj.cc) <- "treatment"

seu.cc.sub.wt <- subset(seu.obj.cc, idents = 'WT')
Idents(seu.cc.sub.wt) <- "seurat_clusters"

seu.cc.sub.ko <- subset(seu.obj.cc, idents = 'GENE-X')
Idents(seu.cc.sub.ko) <- "seurat_clusters"


seu.cc.sub.wt.plot <- DimPlot(seu.cc.sub.wt, label.size = 4.5, label = TRUE) +
  theme(legend.position = "none") + 
  ggtitle("CC-WT-UMAP") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-15,20), ylim = c(-15,20))


seu.cc.sub.ko.plot <- DimPlot(seu.cc.sub.ko, label.size = 4.5, label = TRUE) +
  theme(legend.position = "none") + 
  ggtitle("CC-KO-UMAP") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-15,20), ylim = c(-15,20))


grid.arrange(seu.cc.sub.wt.plot, seu.cc.sub.ko.plot, ncol=2)
```


# Cluster Annotation (CC-Pooled)
```{r}
# annotate CC-Pooled
cluster.id <- c("Microglia", "Microglia", "Endothelial Cells", "Astrocytes",
                "VSMC", "Oligodendrocytes", "Pericytes", "Astrocytes",
                "T Cells", "Mac/Mono", "Astrocytes", "ImmN",
                "Neutrophils", "Schwann Cells", "OPC", "NK Cells",
                "B Cells", "Astrocytes", "VLMC" 
                )

names(cluster.id) <- levels(seu.obj.cc)

seu.obj.cc <- RenameIdents(seu.obj.cc, cluster.id)

# Add 'celltype' column variable
seu.obj.cc@meta.data <- 
  seu.obj.cc@meta.data %>%
  mutate(celltype = seurat_clusters)

seu.obj.cc@meta.data$celltype <- factor(seu.obj.cc@meta.data$celltype, 
                     level = c("0", "1", "2", "3", 
                               "4", "5", "6", "7", 
                               "8", "9", "10", "11", 
                               "12", "13", "14", "15", 
                               "16", "17", "18"),
                     labels = cluster.id)

Idents(seu.obj.cc)  <- "celltype"

seu.obj.cc.plot <- DimPlot(seu.obj.cc, label = TRUE, label.size = 4.5) + 
  theme(legend.position = "none") +
  ggtitle("CC-Pooled-UMAP") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-15,20), ylim = c(-15,20)) +
  theme(axis.line.x = element_line(linewidth = 1.2, color = "black"),
        axis.line.y = element_line(linewidth = 1.2, color = "black"))

feat.plot <- FeaturePlot(seu.obj.cc, features = "Slc1a2") +
  coord_fixed(ratio = 1, xlim = c(-15,20), ylim = c(-15,20))

grid.arrange(seu.obj.cc.plot, feat.plot, ncol=2)
```

# Cluster Annotation (CC-WT-sub)
```{r}
# annotate CC-WT
cluster.id <- c("Microglia", "Microglia", "Endothelial Cells", "Astrocytes", "VSMC", "Oligodendrocytes", "Pericytes", "Astrocytes",
                "T Cells", "Mac/Mono", "Astrocytes", "ImmN",
                "Neutrophils", "Schwann Cells", "OPC", "NK Cells", "B Cells", "Astrocytes", "VLMC")

names(cluster.id) <- levels(seu.cc.sub.wt)

seu.cc.sub.wt <- RenameIdents(seu.cc.sub.wt, cluster.id)

# Add 'celltype' column variable
seu.cc.sub.wt@meta.data <- 
  seu.cc.sub.wt@meta.data %>%
  mutate(celltype = seurat_clusters)

seu.cc.sub.wt@meta.data$celltype <- factor(seu.cc.sub.wt@meta.data$celltype, 
                     level = c("0", "1", "2", "3", 
                               "4", "5", "6", "7", 
                               "8", "9", "10", "11", 
                               "12", "13", "14", "15", 
                               "16", "17", "18"),
                     labels = cluster.id)

Idents(seu.cc.sub.wt)  <- "celltype"
  
seu.cc.sub.wt.plot <- DimPlot(seu.cc.sub.wt, label = TRUE, label.size = 4.5) + 
  theme(legend.position = "none") +
  ggtitle("CC-WT-UMAP") + 
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-15,20), ylim = c(-15,20)) +
  theme(axis.line.x = element_line(linewidth = 1.2, color = "black"),
        axis.line.y = element_line(linewidth = 1.2, color = "black"))

feat.plot <- FeaturePlot(seu.cc.sub.wt, features = "Slc1a2") +
  coord_fixed(ratio = 1, xlim = c(-15,20), ylim = c(-15,20))

grid.arrange(seu.cc.sub.wt.plot, feat.plot, ncol=2)
```


# SAVE Pooled, WT, Seurat object
(for future use by Immunology Team)
```{r, eval = FALSE}

saveRDS(
  seu.obj.cc,
  file = "/hpc/home/jcy13/scRNA_crypto_Jae/Seurat_Objects/seu_obj_cc_020924.rds")

saveRDS(
  seu.cc.sub.wt, 
  file = "/hpc/home/jcy13/scRNA_crypto_Jae/Seurat_Objects/seu_obj_cc_wt_020924.rds")
```

# Violin Plot
```{r}
## remove the x-axis text and tick
## plot.margin to adjust the white space between each plot.
## ... pass any arguments to VlnPlot in Seurat
modify_vlnplot<- function(obj, 
                          feature, 
                          pt.size = 0, 
                          plot.margin = unit(c(-0.75, 0, -0.75, 0), "cm")
                          ) {
  p <- VlnPlot(obj, features = feature, pt.size = pt.size )  + 
  xlab("") + ylab(feature) + ggtitle("") + 
  theme(legend.position = "none",
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.title.y = element_text(size = rel(1), angle = 0),
        axis.text.y = element_text(size = rel(1)),
        plot.margin = plot.margin ) 
  return(p)
}

## extract the max value of the y axis
extract_max<- function(p){
  ymax<- max(ggplot_build(p)$layout$panel_scales_y[[1]]$range$range)
  return(ceiling(ymax))
}

## main function
StackedVlnPlot<- function(obj, features,
                          pt.size = 0, 
                          plot.margin = unit(c(-0.75, 0, -0.75, 0), "cm"),
                          ...) {
  
  plot_list<- purrr::map(features, function(x) modify_vlnplot(obj = obj,feature = x, ...))
  
  # Add back x-axis title to bottom plot. patchwork is going to support this?
  plot_list[[length(plot_list)]]<- plot_list[[length(plot_list)]] +
    theme(axis.text.x=element_text(), axis.ticks.x = element_line())
  
  # change the y-axis tick to only max value 
  ymaxs<- purrr::map_dbl(plot_list, extract_max)
  plot_list<- purrr::map2(plot_list, ymaxs, function(x,y) x + 
                            scale_y_continuous(breaks = c(y)) + 
                            expand_limits(y = y))

  p<- patchwork::wrap_plots(plotlist = plot_list, ncol = 1)
  return(p)
}
```

# Tabulate Cell and Gene counts
## Counting number of *cells* per celltype
```{r}
# (WT)
cell_count_wt <- seu.cc.sub.wt@meta.data %>%
  dplyr::group_by(celltype, time.point) %>%
  dplyr::count(celltype)

cell_count_wt

```

## Counting *genes* per celltype
```{r}
# Wild Type (WT)
# Convert to a dense matrix if not already, then to a data frame
gene_wt_df <- as.data.frame(as.matrix(seu.cc.sub.wt@assays$RNA@counts))

# Transpose the data frame so that genes are rows and cells are columns
gene_wt_t_df <- t(gene_wt_df)
gene_wt_t_df_sparse <- as(gene_wt_t_df, "sparseMatrix")

# Sum unique gene expression counts (i.e. sum number of non-zero instances)
gene_counts_unique <- rowSums(gene_wt_t_df_sparse != 0)

# Sum total number of gene-expression counts per cell
gene_counts_total <- rowSums(gene_wt_t_df_sparse)

# Summarize in df to be inner_joined to metadata df
gene_counts_df <- data.frame(
  unique_counts = gene_counts_unique,
  total_counts = gene_counts_total
  )
row.names(gene_counts_df) <- row.names(gene_wt_t_df_sparse)
gene_counts_df <- tibble::rownames_to_column(gene_counts_df, var = "cell_barcode")
rm(gene_wt_df, gene_wt_t_df_sparse, gene_counts_unique, gene_counts_total)  

# inner join gene_counts_df to seu.cc.sub.wt metadata
meta_gene_wt_df <- seu.cc.sub.wt@meta.data %>%
  tibble::rownames_to_column(var="cell_barcode") %>%
  dplyr::inner_join(., gene_counts_df, by = "cell_barcode")

gene_unique_count_wt <- meta_gene_wt_df %>%
  dplyr::group_by(celltype, time.point) %>%
  dplyr::count(celltype, wt = unique_counts)

gene_total_count_wt <- meta_gene_wt_df %>%
  dplyr::group_by(celltype, time.point) %>%
  dplyr::count(celltype, wt = total_counts)

gene_unique_count_wt
gene_total_count_wt
```

# [Request #3]
```{r}
gene_list <- c("Tmem119", "P2ry12", "Gpr34", "Sall1", 
               "Aqp4", "Gfap", "Gja1", "Slc1a2", "Aldh1l1",
               "Cldn5", "Rax",
               "Acta2", "Cspg4",
               "Vtn", "Kcnj8",
               "Pdgfra",
               "Cd3e",
               "Lyz2", "Plac8", "Cd86", "Adgre1", "Cd40",
               "Sox11", "Dcx",
               "S100a9", "Ly6g",
               "Npy",
               "Igkc", "Cd79a", "Cd19",
               "Slc6a13",
               "Folr1",
               "Nkg7"
               )


Idents(seu.cc.sub.wt) <- "celltype"

dot.plot.1 <- DotPlot(seu.cc.sub.wt, features = gene_list) +
  coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust=1))
```

## Dot Plot
```{r}
gene_list <- c("Tmem119", "P2ry12", "Gpr34", "Sall1", 
               "Aqp4", "Gfap", "Gja1", "Slc1a2", "Aldh1l1",
               "Cldn5", "Rax",
               "Acta2", "Cspg4",
               "Vtn", "Kcnj8",
               "Pdgfra",
               "Cd3e",
               "Lyz2", "Plac8", "Cd86", "Adgre1", "Cd40",
               "Sox11", "Dcx",
               "S100a9", "Ly6g",
               "Npy",
               "Igkc", "Cd79a", "Cd19",
               "Slc6a13",
               "Folr1",
               "Nkg7"
               )

# new celltype ordering
new.lvl.order <- c(
  "Microglia",
  "Astrocytes",
  "Endothelial Cells",
  "VSMC",
  "Pericytes",
  "Oligodendrocytes",
  "T Cells",
  "Mac/Mono",
  "OPC",
  "ImmN",
  "Neutrophils",
  "Schwann Cells",
  "B Cells",
  "VLMC",
  "NK Cells"
  )

# create seurat clone object v2
seu.cc.sub.wt.2 <- seu.cc.wt

# assign new level ordering to seu obj v2
seu.cc.sub.wt.2@meta.data$celltype <- factor(seu.cc.sub.wt.2@meta.data$celltype,
                                             levels = new.lvl.order
                                           )

# dot-plot plotting
Idents(seu.cc.sub.wt.2) <- "celltype"
dot.plot.wt <- DotPlot(seu.cc.sub.wt.2, features = gene_list) +
  coord_flip() + 
  theme(axis.text.x = element_text(angle = 45, hjust=1))

rm(seu.cc.sub.wt.2)
dot.plot.wt
```

# [Disregard] Horizontal Bar Chart (seu.obj.cc.wt)
```{r, eval = FALSE}
# Summary data

# (A) get total counts by cell type in order of levels defined for 'celltype'
cell_tot <- as.data.frame(table(seu.cc.sub.wt@meta.data$celltype))

# (B) get total count of genes detected per cell type
# extract the expression matrix from the Seurat object
expr_matrix <- seu.cc.sub.wt[["RNA"]]@counts 

# If you're working with a sparse matrix, convert it to a dense matrix (this will be necessary for the next steps)
if (class(expr_matrix) == "dgCMatrix") {
  expr_matrix <- as.matrix(expr_matrix)
}

# Initialize a vector to store the number of genes per cell type
genes_per_celltype <- numeric()

# Get the list of cell types
cell_types <- levels(seu.cc.sub.wt@meta.data$celltype)

# Calculate number of genes for each cell type
for (ct in cell_types) {
  Cells_of_ct <- which(seu.cc.sub.wt@meta.data$celltype == ct)
  subset_matrix <- expr_matrix[, Cells_of_ct]
  detected_genes <- sum(apply(subset_matrix, 1, function(row) any(row > 0)))
  genes_per_celltype[ct] <- detected_genes
}

gene_df <- as.data.frame(genes_per_celltype)

# (C) append all information above
data <- data.frame(
  celltype = levels(seu.cc.sub.wt@meta.data$celltype),
  cell_count = cell_tot[,2],
  gene_count = gene_df[,1]
  # type = rep(c("Cell Count", "Gene Count"), times = 6)
)

long.data <- tidyr::pivot_longer(data, cols = cell_count:gene_count)

# Plot for Cell Count
p1 <- ggplot(subset(long.data, name == "cell_count")) +
  geom_text(aes(x = celltype, y = 2, label = value), hjust = 0.25) +
  # Adjust expansion here
  scale_y_discrete(expand = expansion(mult = c(0.05, 0.05))) + 
  coord_flip() +
  labs(y = "No. of Cells", x = "") +
  theme_minimal() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.x = element_blank())

# Plot for Gene Count
p2 <- ggplot(subset(long.data, name == "gene_count")) +
  geom_bar(aes(x = celltype, y = value, fill = celltype), stat = "identity", position = "dodge", width = 0.9) +
  coord_flip() +
  labs(y = "Detected genes per cell type", x = "") +
  theme_minimal() +
  scale_fill_brewer(palette = "Set3") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        axis.text.y = element_blank(),
        legend.position = "none")

# Combine plots
grid.arrange(p1, p2, ncol=2)
```


# [Request #2]

## Microglia (WT)
```{r}
mg.feat.list <- c("P2ry12", "Gpr34", "Tmem119", "Sall1")

mg.plot <- list()

# iterate each gene plot in order to control final layout (legend display, axis, etc.)
for (i in 1:length(mg.feat.list)){
  plot <- FeaturePlot(seu.cc.sub.wt, 
                      features = mg.feat.list[i], 
                      keep.scale = "all",
                      cols = c("lightgrey", "#f02602")
                      ) + 
    theme(legend.text = element_blank(),
          legend.key.height = unit(0.4, "cm"),   
          legend.key.width = unit(0.4, "cm")
          ) +    
    coord_fixed(ratio = 1, xlim = c(-15,20), ylim = c(-15,20))
    
  # if the plot is not the last one, remove the legend scale
  if(i != length(mg.feat.list)){
    plot <- plot + theme(legend.position = "none", legend.text = element_blank())
  }
  
  # final gene plot list
  mg.plot[[i]] <- plot
}

# display gene plots with patchwork-plot layout function
combine.mg.plot <- mg.plot[[1]] + mg.plot[[2]] + mg.plot[[3]] + mg.plot[[4]] +
  patchwork::plot_layout(ncol=4, nrow=1)

combine.mg.plot
```

## Microglia (WT) [ISO]
```{r}
# subset out Microglia
Idents(seu.cc.sub.wt) <- "celltype"
wt.mg <- subset(seu.cc.sub.wt, idents="Microglia")


mg.feat.list <- c("P2ry12", "Gpr34", "Tmem119", "Sall1")
num_col <- length(mg.feat.list)

mg.plot <- list()

# iterate each gene plot in order to control final layout (legend display, axis, etc.)
for (i in 1:length(mg.feat.list)){
  plot <- FeaturePlot(wt.mg, 
                      features = mg.feat.list[i], 
                      keep.scale = "all",
                      cols = c("lightgrey", "#f02602"),
                      pt.size = 0.5
                      ) + 
    theme(legend.text = element_blank(),
          legend.key.height = unit(0.4, "cm"),   
          legend.key.width = unit(0.4, "cm")
          ) +    
    coord_fixed(ratio = 1, xlim = c(-15,5), ylim = c(-15,5))
    
  # if the plot is not the last one, remove the legend scale
  if(i != length(mg.feat.list)){
    plot <- plot + theme(legend.position = "none", legend.text = element_blank())
  }
  
  # final gene plot list
  mg.plot[[i]] <- plot
}

# display gene plots with patchwork-plot layout function
combine.mg.plot <- mg.plot[[1]] + mg.plot[[2]] + mg.plot[[3]] + mg.plot[[4]] +
  patchwork::plot_layout(ncol=num_col, nrow=1)

combine.mg.plot
```

## Astrocytes (WT)
```{r}
astro.feat.list <- c("Slc1a2", "Aqp4", "Gfap", "Adgrv1")

astro.list <- list()

# iterate each gene plot in order to control final layout (legend display, axis, etc.)
for (i in 1:length(astro.feat.list)){
  plot <- FeaturePlot(seu.cc.sub.wt, 
                      features = astro.feat.list[i], 
                      keep.scale = "all",
                      cols = c("lightgrey", "#f02602")
                      ) + 
    theme(legend.text = element_blank(),
          legend.key.height = unit(0.4, "cm"),   
          legend.key.width = unit(0.4, "cm")
          ) +    
    coord_fixed(ratio = 1, xlim = c(-15,20), ylim = c(-15,20))
    
  # if the plot is not the last one, remove the legend scale
  if(i != length(astro.feat.list)){
    plot <- plot + theme(legend.position = "none", legend.text = element_blank())
  }
  
  # final gene plot list
  astro.list[[i]] <- plot
}

# display gene plots with patchwork-plot layout function
combine.astro.plot <- astro.list[[1]] + astro.list[[2]] + astro.list[[3]] + astro.list[[4]] +
  patchwork::plot_layout(ncol=4, nrow=1)

```

## Astrocytes (WT) [ISO]
```{r}
# subset out Astrocytes
Idents(seu.cc.sub.wt) <- "celltype"
wt.astro <- subset(seu.cc.sub.wt, idents="Astrocytes")


astro.feat.list <- c("Slc1a2", "Aqp4", "Gfap", "Adgrv1")

astro.list <- list()

# iterate each gene plot in order to control final layout (legend display, axis, etc.)
for (i in 1:length(astro.feat.list)){
  plot <- FeaturePlot(wt.astro, 
                      features = astro.feat.list[i], 
                      keep.scale = "all",
                      cols = c("lightgrey", "#f02602"),
                      pt.size = 0.5
                      ) + 
    theme(legend.text = element_blank(),
          legend.key.height = unit(0.4, "cm"),   
          legend.key.width = unit(0.4, "cm")
          ) +    
    coord_fixed(ratio = 1, xlim = c(0,20), ylim = c(-15,5))
    
  # if the plot is not the last one, remove the legend scale
  if(i != length(astro.feat.list)){
    plot <- plot + theme(legend.position = "none", legend.text = element_blank())
  }
  
  # final gene plot list
  astro.list[[i]] <- plot
}

# display gene plots with patchwork-plot layout function
combine.astro.plot <- astro.list[[1]] + astro.list[[2]] + astro.list[[3]] + astro.list[[4]] +
  patchwork::plot_layout(ncol=4, nrow=1)

```

## Oligodendrocytes (WT)
```{r}
oligo.feat.list <- c("Cldn11", "Pdgfra", "Sox10", "Olig1")

oligo.plot <- list()

# iterate each gene plot in order to control final layout (legend display, axis, etc.)
for (i in 1:length(oligo.feat.list)){
  plot <- FeaturePlot(seu.cc.sub.wt, 
                      features = oligo.feat.list[i], 
                      keep.scale = "all",
                      cols = c("lightgrey", "#f02602")
                      ) + 
    theme(legend.text = element_blank(),
          legend.key.height = unit(0.4, "cm"),   
          legend.key.width = unit(0.4, "cm")
          ) +    
    coord_fixed(ratio = 1, xlim = c(-15,20), ylim = c(-15,20))
    
  # if the plot is not the last one, remove the legend scale
  if(i != length(oligo.feat.list)){
    plot <- plot + theme(legend.position = "none", legend.text = element_blank())
  }
  
  # final gene plot list
  oligo.plot[[i]] <- plot
}

# display gene plots with patchwork-plot layout function
combine.oligo.plot <- oligo.plot[[1]] + oligo.plot[[2]] + oligo.plot[[3]] + oligo.plot[[4]] +
  patchwork::plot_layout(ncol=4, nrow=1)

```

## Oligodendrocytes (WT) [ISO]
```{r}
# subset out Oligo
Idents(seu.cc.sub.wt) <- "celltype"
wt.oligo <- subset(seu.cc.sub.wt, idents="Oligodendrocytes")

oligo.feat.list <- c("Cldn11", "Pdgfra", "Sox10", "Olig1")
oligo.plot <- list()

# iterate each gene plot in order to control final layout (legend display, axis, etc.)
for (i in 1:length(oligo.feat.list)){
  plot <- FeaturePlot(wt.oligo, 
                      features = oligo.feat.list[i], 
                      keep.scale = "all",
                      cols = c("lightgrey", "#f02602"),
                      pt.size = 0.5
                      ) + 
    theme(legend.text = element_blank(),
          legend.key.height = unit(0.4, "cm"),   
          legend.key.width = unit(0.4, "cm")
          ) +    
    coord_fixed(ratio = 1, xlim = c(0,15), ylim = c(-20,-5))
    
  # if the plot is not the last one, remove the legend scale
  if(i != length(oligo.feat.list)){
    plot <- plot + theme(legend.position = "none", legend.text = element_blank())
  }
  
  # final gene plot list
  oligo.plot[[i]] <- plot
}

# display gene plots with patchwork-plot layout function
combine.oligo.plot <- oligo.plot[[1]] + oligo.plot[[2]] + oligo.plot[[3]] + oligo.plot[[4]] +
  patchwork::plot_layout(ncol=4, nrow=1)

```

# Celltype UMAP
## Microglia
```{r}
Idents(seu.obj.cc) <- "celltype"
seu.obj.cc.mg <- subset(seu.obj.cc, idents="Microglia")

seu.obj.cc.mg <- FindVariableFeatures(seu.obj.cc.mg, selection.method = "vst")
seu.obj.cc.mg <- ScaleData(seu.obj.cc.mg, selection.method = "vst")
seu.obj.cc.mg <- RunPCA(seu.obj.cc.mg, verbose = FALSE)
#ElbowPlot(seu.obj.cc.mg, ndims = 50, reduction = "pca")

# compute dimension of k.param by taking square root of cell sample count
n = sqrt(seu.obj.cc.mg@assays$RNA@counts@Dim[2])
seu.obj.cc.mg <- Seurat::FindNeighbors(seu.obj.cc.mg, 
                                    k.param = n, 
                                    dims = 1:20
                                    )

seu.obj.cc.mg <- Seurat::FindClusters(seu.obj.cc.mg,
                                   algorithm = 1,
                                   resolution = 0.4,
                                   random.seed = 72423
                                   )

seu.obj.cc.mg <- RunUMAP(seu.obj.cc.mg, reduction = "pca", dims = 1:20)

Seurat::DefaultAssay(seu.obj.cc.mg) <- "RNA"
Idents(seu.obj.cc.mg) <- "seurat_clusters"


seu.obj.cc.mg.plot <- DimPlot(seu.obj.cc.mg, label = TRUE, label.size = 3) + 
  ggtitle("CC-Microglia-Pooled-UMAP") +
  theme(plot.title = element_text(hjust = 0.5)) + 
  coord_fixed(ratio = 1, xlim = c(-15, 10), ylim = c(-15, 10))

feat.plot <- FeaturePlot(seu.obj.cc.mg, features = "Cx3cr1") +
  coord_fixed(ratio = 1, xlim = c(-15, 10), ylim = c(-15, 10))

grid.arrange(seu.obj.cc.mg.plot, feat.plot, ncol=2)
Idents(seu.obj.cc.mg) <- "seurat_clusters"
```

```{r}
Idents(seu.obj.cc.mg) <- "treatment"
seu.obj.cc.mg.wt <- subset(seu.obj.cc.mg, idents="WT")

# --------------------------------------------------------------------------------------------------------
# Note: following lines of code (i.e. re-ordering cluster labels cosmetically impact pseudotime trajectory plots)
# re-order labels for clusters (e.g. 0, 1, 2) for more intuitive pseudo-time trajectory plot interpretation further down below 
# 2 -> 0
# 0 -> 1
# 1 -> 2
seu.obj.cc.mg.wt@meta.data <- seu.obj.cc.mg.wt@meta.data %>%
  dplyr::mutate(seurat_clusters = as.numeric(as.character(seurat_clusters))) %>%
  dplyr::mutate(seurat_clusters = seurat_clusters + 1) %>%
  dplyr::mutate(seurat_clusters = ifelse(seurat_clusters == 3, 0, seurat_clusters))

# --------------------------------------------------------------------------------------------------------

Idents(seu.obj.cc.mg.wt) <- "seurat_clusters"
seu.obj.cc.mg.wt.plot <- DimPlot(seu.obj.cc.mg.wt, label = TRUE, label.size = 3) + 
  ggtitle("CC-Microglia-WT-UMAP") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-15,10), ylim = c(-15,10))

seu.obj.cc.mg.wt.plot
```

### Pseudotime Trajectory (WT) [Request #6]
*important notes* 
- currently do not have labels or identifiers for the 3x cell clusters identified in Microglia
- when determining which cluster group (i.e. 0, 1, 2) is root, trajectory was first ran with left group (#1) set to root. From time panel plotting we know that left group #1 is actually the final group because it is all cells in day 14 of time panel plot. Thus, we can deduce that whichever cluster group is identified as the final cluster (#<insert>) from this run is actually the true root cluster group since trajectory graphs are *unidirectional*. 

(i) convert to cell_data_set object
```{r}
# Online recommended "as.CellDataSet()" function DOES NOT seem to successfully convert 
# Seurat to Monocle3 CDS object
# *Manually* create *new* CDS object from scratch (i.e. original matrices) based on link:
# https://github.com/satijalab/seurat/issues/2833

# get expression matrix
temp.mat <- seu.obj.cc.mg.wt@assays[["RNA"]]@counts
temp.mat <- temp.mat[rownames(seu.obj.cc.mg.wt@reductions[["pca"]]@feature.loadings), ]
express.mat <- temp.mat

# get cell_metadata matrix
cell.meta <- as.data.frame(
  seu.obj.cc.mg.wt@assays[["RNA"]]@counts@Dimnames[[2]],
  row.names = seu.obj.cc.mg.wt@assays[["RNA"]]@counts@Dimnames[[2]]
  )
colnames(cell.meta) <- "barcode"

# get gene_metadata matrix
gene.mat <- as.data.frame(
  rownames(seu.obj.cc.mg.wt@reductions[["pca"]]@feature.loadings),
  row.names = rownames(seu.obj.cc.mg.wt@reductions[["pca"]]@feature.loadings)
  )
colnames(gene.mat) <- "gene_short_name"

# declare *new* CDS object
mg.wt.cds <- new_cell_data_set(
  express.mat, 
  cell_metadata = cell.meta, 
  gene_metadata = gene.mat
  )

recreate.partition <- c(rep(1, length(mg.wt.cds@colData@rownames)))
names(recreate.partition) <- mg.wt.cds@colData@rownames
recreate.partition <- as.factor(recreate.partition)

mg.wt.cds@clusters@listData[["UMAP"]][["partitions"]] <- recreate.partition
```

(ii) cluster cells (based on seurat UMAP)
```{r}
# map over additional seurat info to new CDS object
list_cluster <- seu.obj.cc.mg.wt@active.ident
names(list_cluster) <- seu.obj.cc.mg.wt@assays[["RNA"]]@data@Dimnames[[2]]

mg.wt.cds@clusters@listData[["UMAP"]][["clusters"]] <- list_cluster
mg.wt.cds@clusters@listData[["UMAP"]][["louvain_res"]] <- "NA"
mg.wt.cds@int_colData@listData$reducedDims@listData[["UMAP"]] <- seu.obj.cc.mg.wt@reductions[["umap"]]@cell.embeddings
```

(iii) learn trajectory graph
```{r}
mg.wt.cds <- learn_graph(mg.wt.cds, use_partition = FALSE)

ymin = -15
ymax = 10

pre_order.plot <- plot_cells(
  mg.wt.cds,
  color_cells_by = 'cluster',
  cell_size = 0.6,
  label_groups_by_cluster = TRUE,
  group_label_size = 5,
  label_branch_points = FALSE,
  label_roots = FALSE
) +
  theme(legend.position = "right") +
  scale_y_continuous(limits = c(ymin, ymax))
```

(iv) order cells by pseudotime
```{r}
# order_cells() was first ran with group #2 set as root
# from our time plots we know that group #2 is actually the final pseudotime group
# so we can infer that the true root group is #0
# and final estimated trajectory is #0 —> #1 —> #2
mg.wt.cds <- order_cells(
  mg.wt.cds,
  reduction_method = 'UMAP',
  root_cells = colnames(mg.wt.cds[,clusters(mg.wt.cds) == '0'])
  )

ymin = -15
ymax = 10

post_order.plot <- plot_cells(
  mg.wt.cds,
  color_cells_by = 'pseudotime',
  cell_size = 0.6,
  group_label_size = 10,
  label_groups_by_cluster = FALSE,
  label_branch_points = FALSE,
  label_leaves = FALSE,
  label_roots = FALSE
) +
  theme(legend.position = "right") +
  scale_y_continuous(limits = c(ymin, ymax))

# # boxplot representation of pseudotime trajectory
# mg.wt.cds$mono3_pseudotime <- pseudotime(mg.wt.cds)
# data.pseudo <- as.data.frame(pData(mg.wt.cds))
# 
# # run boxplot through ggplot
# ggplot(data.pseudo, aes(mono3_pseudotime, reorder(cluster, mono3_pseudotime, median), fill = cluster)) +
#   geom_boxplot()
```

(v) visualize pseudotime with seurat
```{r}
ymin = -15
ymax = 10

xmin = -15
xmax = 10

# visualizing pseudotime in seurat
seu.obj.cc.mg.wt$pseudotime <- pseudotime(mg.wt.cds)
Idents(seu.obj.cc.mg.wt) <- seu.obj.cc.mg.wt$seurat_clusters
pseudo_feature_plot <- FeaturePlot(
  seu.obj.cc.mg.wt, features = "pseudotime", label = T, pt.size = 0.6) + 
  scale_y_continuous(limits = c(ymin, ymax)) +
  scale_x_continuous(limits = c(xmin, xmax))
```

(vi) genes that change as function of pseudotime
```{r}
# Finding genes that change as a function of pseudotime
deg.mg.cells <- graph_test(
  mg.wt.cds, 
  neighbor_graph = 'principal_graph', 
  cores = 4
  )

deg.mg.cells %>% 
  arrange(q_value) %>% 
  filter(status == 'OK') %>% 
  head()

FeaturePlot(seu.obj.cc.mg.wt, features = c('Cxcl9', 'Ccl5', 'Fabp5'))
```

### Cell Cycle Plot (WT) [Request #5] 
```{r}
# Order factor levels (G1, S, G2M) to align with desired colors
phase_colors <- brewer.pal(3, "Set1")
seu.obj.cc.mg.wt@meta.data$Phase <- factor(
  seu.obj.cc.mg.wt@meta.data$Phase, 
  levels = c("S", "G1", "G2M")
  )
all_phases <- levels(seu.obj.cc.mg.wt@meta.data$Phase)

# Check for length misalignment
if(length(all_phases) > length(phase_colors)) {
  stop("Number of phases exceeds length of color palette. Choose larger palette or adjust accordingly.")
}
# Fix colors for cell-cycle phases (G1=blue, S=red, G2M=green)
names(phase_colors) <- all_phases

Idents(seu.obj.cc.mg.wt) <- "Phase"
seu.obj.cc.mg.wt.plot <- DimPlot(
  seu.obj.cc.mg.wt,
  label = FALSE, 
  label.size = 3,
  cols = alpha(phase_colors, 0.6),
  group.by = "Phase"
  ) + 
  ggtitle("CC-Microglia-WT-UMAP") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-15,10), ylim = c(-15,10))

seu.obj.cc.mg.wt.plot
```

## Astrocytes
```{r}
Idents(seu.obj.cc) <- "celltype"
seu.obj.cc.astro <- subset(seu.obj.cc, idents="Astrocytes")

seu.obj.cc.astro <- FindVariableFeatures(seu.obj.cc.astro, selection.method = "vst")
seu.obj.cc.astro <- ScaleData(seu.obj.cc.astro, selection.method = "vst")
seu.obj.cc.astro <- RunPCA(seu.obj.cc.astro, verbose = FALSE)
ElbowPlot(seu.obj.cc.astro, ndims = 50, reduction = "pca")

n = sqrt(seu.obj.cc.astro@assays$RNA@counts@Dim[2])
seu.obj.cc.astro <- Seurat::FindNeighbors(seu.obj.cc.astro, 
                                    k.param = n, 
                                    dims = 1:20
                                    )

seu.obj.cc.astro <- Seurat::FindClusters(seu.obj.cc.astro,
                                   algorithm = 1,
                                   resolution = 0.1,
                                   random.seed = 72423
                                   )

seu.obj.cc.astro <- RunUMAP(seu.obj.cc.astro, reduction = "pca", dims = 1:30)

Seurat::DefaultAssay(seu.obj.cc.astro) <- "RNA"
Idents(seu.obj.cc.astro) <- "seurat_clusters"


seu.obj.cc.astro.plot <- DimPlot(seu.obj.cc.astro, label = TRUE, label.size = 3) + 
  ggtitle("CC-Astrocytes-Pooled-UMAP") +
  theme(plot.title = element_text(hjust = 0.5)) + 
  coord_fixed(ratio = 1, xlim = c(-10,30), ylim = c(-20,20))

feat.plot <- FeaturePlot(seu.obj.cc.astro, features = "P2ry12")

grid.arrange(seu.obj.cc.astro.plot, feat.plot, ncol=2)
```


```{r}
Idents(seu.obj.cc.astro) <- "treatment"
seu.obj.cc.astro.wt <- subset(seu.obj.cc.astro, idents="WT")

# --------------------------------------------------------------------------------------------------------
# Note: following lines of code (i.e. re-ordering cluster labels cosmetically impact pseudotime trajectory plots)
# re-order labels for clusters (e.g. 0, 1, 2) for more intuitive pseudo-time trajectory plot interpretation further down below 
# 0 -> 0
# 1 -> 5
# 2 -> 4
# 3 -> 3
# 4 -> 2
# 5 -> 1
# 6 -> 6

seu.obj.cc.astro.wt@meta.data <- seu.obj.cc.astro.wt@meta.data %>%
  dplyr::mutate(seurat_clusters = as.numeric(as.character(seurat_clusters))) %>%
  dplyr::mutate(seurat_clusters = case_when(
    seurat_clusters == 0 ~ 0,
    seurat_clusters == 1 ~ 5,
    seurat_clusters == 2 ~ 4,
    seurat_clusters == 3 ~ 3,
    seurat_clusters == 4 ~ 2,
    seurat_clusters == 5 ~ 1,
    seurat_clusters == 6 ~ 6
  ))
# --------------------------------------------------------------------------------------------------------


Idents(seu.obj.cc.astro.wt) <- "seurat_clusters"
seu.obj.cc.astro.wt.plot <- DimPlot(seu.obj.cc.astro.wt, label = TRUE, label.size = 3) + 
  ggtitle("CC-Astrocytes-WT-UMAP") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,30), ylim = c(-20,20))

seu.obj.cc.astro.wt.plot
```

### Pseudotime Trajectory (WT) [Request #6]
```{r}
# seu.obj.cc.astro.wt
ymin = -20
ymax = 20
xmin = -10
xmax = 30


# (i) convert to cell_data_set object
# Online recommended "as.CellDataSet()" function DOES NOT seem to successfully convert 
# Seurat to Monocle3 CDS object
# *Manually* create *new* CDS object from scratch (i.e. original matrices) based on link:
# https://github.com/satijalab/seurat/issues/2833

# get expression matrix
temp.mat <- seu.obj.cc.astro.wt@assays[["RNA"]]@counts
temp.mat <- temp.mat[rownames(seu.obj.cc.astro.wt@reductions[["pca"]]@feature.loadings), ]
express.mat <- temp.mat

# get cell_metadata matrix
cell.meta <- as.data.frame(
  seu.obj.cc.astro.wt@assays[["RNA"]]@counts@Dimnames[[2]],
  row.names = seu.obj.cc.astro.wt@assays[["RNA"]]@counts@Dimnames[[2]]
  )
colnames(cell.meta) <- "barcode"

# get gene_metadata matrix
gene.mat <- as.data.frame(
  rownames(seu.obj.cc.astro.wt@reductions[["pca"]]@feature.loadings),
  row.names = rownames(seu.obj.cc.astro.wt@reductions[["pca"]]@feature.loadings)
  )
colnames(gene.mat) <- "gene_short_name"

# declare *new* CDS object
astro.wt.cds <- new_cell_data_set(
  express.mat, 
  cell_metadata = cell.meta, 
  gene_metadata = gene.mat
  )

recreate.partition <- c(rep(1, length(astro.wt.cds@colData@rownames)))
names(recreate.partition) <- astro.wt.cds@colData@rownames
recreate.partition <- as.factor(recreate.partition)

astro.wt.cds@clusters@listData[["UMAP"]][["partitions"]] <- recreate.partition


# (ii) cluster cells (based on seurat UMAP)
# map over additional seurat info to new CDS object
list_cluster <- seu.obj.cc.astro.wt@active.ident
names(list_cluster) <- seu.obj.cc.astro.wt@assays[["RNA"]]@data@Dimnames[[2]]

astro.wt.cds@clusters@listData[["UMAP"]][["clusters"]] <- list_cluster
astro.wt.cds@clusters@listData[["UMAP"]][["louvain_res"]] <- "NA"
astro.wt.cds@int_colData@listData$reducedDims@listData[["UMAP"]] <- seu.obj.cc.astro.wt@reductions[["umap"]]@cell.embeddings


# (iii) learn trajectory graph
astro.wt.cds <- learn_graph(astro.wt.cds, use_partition = FALSE)

pre_order.plot <- plot_cells(
  astro.wt.cds,
  color_cells_by = 'cluster',
  cell_size = 0.6,
  label_groups_by_cluster = TRUE,
  group_label_size = 5,
  label_branch_points = FALSE,
  label_roots = FALSE
) +
  theme(legend.position = "right") +
  scale_y_continuous(limits = c(ymin, ymax))


# (iv) order cells by pseudotime
astro.wt.cds <- order_cells(
  astro.wt.cds,
  reduction_method = 'UMAP',
  root_cells = colnames(astro.wt.cds[,clusters(astro.wt.cds) == '0'])
  )

post_order.plot <- plot_cells(
  astro.wt.cds,
  color_cells_by = 'pseudotime',
  cell_size = 0.6,
  group_label_size = 10,
  label_groups_by_cluster = FALSE,
  label_branch_points = FALSE,
  label_leaves = FALSE,
  label_roots = FALSE
) +
  theme(legend.position = "right") +
  scale_y_continuous(limits = c(ymin, ymax))

# (v) genes that change as function of pseudotime
# visualizing pseudotime in seurat
seu.obj.cc.astro.wt$pseudotime <- pseudotime(astro.wt.cds)
Idents(seu.obj.cc.astro.wt) <- seu.obj.cc.astro.wt$seurat_clusters
pseudo_feature_plot <- FeaturePlot(
  seu.obj.cc.astro.wt, features = "pseudotime", label = T, pt.size = 0.6) +
  scale_y_continuous(limits = c(ymin, ymax)) +
  scale_x_continuous(limits = c(xmin, xmax))
```

### Cell Cycle Plot (WT) [Request #5]
```{r}
# Order factor levels (G1, S, G2M) to align with desired colors
phase_colors <- brewer.pal(3, "Set1")
seu.obj.cc.astro.wt@meta.data$Phase <- factor(
  seu.obj.cc.astro.wt@meta.data$Phase,
  levels = c("S", "G1", "G2M")
)
all_phases <- levels(seu.obj.cc.astro.wt@meta.data$Phase)

# Check for length misalignment
if(length(all_phases) > length(phase_colors)){
  stop("Number of phases exceeds length of color palette. Choose larger palette or adjust accordingly.")
}

# Fix colors for cell-cycle phases (G1=blue, S=red, G2M=green)
names(phase_colors) <- all_phases

Idents(seu.obj.cc.astro.wt) <- "Phase"
seu.obj.cc.astro.wt.plot <- DimPlot(
  seu.obj.cc.astro.wt, 
  label = FALSE, 
  label.size = 3,
  cols = alpha(phase_colors, 0.6),
  group.by = "Phase"
  ) + 
  ggtitle("CC-Astrocytes-WT-UMAP") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-20,10))

seu.obj.cc.astro.wt.plot
```

## Oligodendrocytes
```{r}
Idents(seu.obj.cc) <- "celltype"
seu.obj.cc.oligo <- subset(seu.obj.cc, idents="Oligodendrocytes")

seu.obj.cc.oligo <- FindVariableFeatures(seu.obj.cc.oligo, selection.method = "vst")
seu.obj.cc.oligo <- ScaleData(seu.obj.cc.oligo, selection.method = "vst")
seu.obj.cc.oligo <- RunPCA(seu.obj.cc.oligo, verbose = FALSE)
ElbowPlot(seu.obj.cc.oligo, ndims = 50, reduction = "pca")

n = sqrt(seu.obj.cc.oligo@assays$RNA@counts@Dim[2])
seu.obj.cc.oligo <- Seurat::FindNeighbors(seu.obj.cc.oligo, 
                                    k.param = n, 
                                    dims = 1:20
                                    )

seu.obj.cc.oligo <- Seurat::FindClusters(seu.obj.cc.oligo,
                                   algorithm = 1,
                                   resolution = 0.3,
                                   random.seed = 72423
                                   )

seu.obj.cc.oligo <- RunUMAP(seu.obj.cc.oligo, reduction = "pca", dims = 1:30)

Seurat::DefaultAssay(seu.obj.cc.oligo) <- "RNA"
Idents(seu.obj.cc.oligo) <- "seurat_clusters"


seu.obj.cc.oligo.plot <- DimPlot(seu.obj.cc.oligo, label = TRUE, label.size = 3) + 
  ggtitle("CC-Oligodendrocytes-Pooled-UMAP") +
  theme(plot.title = element_text(hjust = 0.5)) + 
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-15,15))

feat.plot <- FeaturePlot(seu.obj.cc.oligo, features = "P2ry12")

grid.arrange(seu.obj.cc.oligo.plot, feat.plot, ncol=2)

# cc.mg.markers.1 <- FindMarkers(seu.obj.cc.oligo, ident.1 = 1)
# cc.mg.markers.2 <- FindMarkers(seu.obj.cc.oligo, ident.1 = 2)
# cc.mg.markers.3 <- FindMarkers(seu.obj.cc.oligo, ident.1 = 3)
# cc.mg.markers.4 <- FindMarkers(seu.obj.cc.oligo, ident.1 = 4)
# cc.mg.markers.5 <- FindMarkers(seu.obj.cc.oligo, ident.1 = 5)
# cc.mg.markers.6 <- FindMarkers(seu.obj.cc.oligo, ident.1 = 6)
# cc.mg.markers.7 <- FindMarkers(seu.obj.cc.oligo, ident.1 = 7)
```



```{r}
Idents(seu.obj.cc.oligo) <- "treatment"
seu.obj.cc.oligo.wt <- subset(seu.obj.cc.oligo, idents="WT")

# --------------------------------------------------------------------------------------------------------
# Note: following lines of code (i.e. re-ordering cluster labels cosmetically impact pseudotime trajectory plots)
# re-order labels for clusters (e.g. 0, 1, 2) for more intuitive pseudo-time trajectory plot interpretation further down below 
# 0 -> 0
# 4 -> 1
# 3 -> 2
# 2 -> 3
# 1 -> 4
seu.obj.cc.oligo.wt@meta.data <- seu.obj.cc.oligo.wt@meta.data %>%
  dplyr::mutate(seurat_clusters = as.numeric(as.character(seurat_clusters))) %>%
  dplyr::mutate(seurat_clusters = case_when(
    seurat_clusters == 0 ~ 0,
    seurat_clusters == 1 ~ 4,
    seurat_clusters == 2 ~ 3,
    seurat_clusters == 3 ~ 2,
    seurat_clusters == 4 ~ 1
  ))

# --------------------------------------------------------------------------------------------------------

Idents(seu.obj.cc.oligo.wt) <- "seurat_clusters"
seu.obj.cc.oligo.wt.plot <- DimPlot(seu.obj.cc.oligo.wt, label = TRUE, label.size = 3) + 
  ggtitle("CC-Oligodendrocytes-WT-UMAP") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-15,15))

seu.obj.cc.oligo.wt.plot
```

### Pseudotime Trajectory (WT) [Request #6]
```{r}
# seu.obj.cc.oligo.wt
ymin = -20
ymax = 10
xmin = -10
xmax = 20


# (i) convert to cell_data_set object
# Online recommended "as.CellDataSet()" function DOES NOT seem to successfully convert 
# Seurat to Monocle3 CDS object
# *Manually* create *new* CDS object from scratch (i.e. original matrices) based on link:
# https://github.com/satijalab/seurat/issues/2833

# get expression matrix
temp.mat <- seu.obj.cc.oligo.wt@assays[["RNA"]]@counts
temp.mat <- temp.mat[rownames(seu.obj.cc.oligo.wt@reductions[["pca"]]@feature.loadings), ]
express.mat <- temp.mat

# get cell_metadata matrix
cell.meta <- as.data.frame(
  seu.obj.cc.oligo.wt@assays[["RNA"]]@counts@Dimnames[[2]],
  row.names = seu.obj.cc.oligo.wt@assays[["RNA"]]@counts@Dimnames[[2]]
  )
colnames(cell.meta) <- "barcode"

# get gene_metadata matrix
gene.mat <- as.data.frame(
  rownames(seu.obj.cc.oligo.wt@reductions[["pca"]]@feature.loadings),
  row.names = rownames(seu.obj.cc.oligo.wt@reductions[["pca"]]@feature.loadings)
  )
colnames(gene.mat) <- "gene_short_name"

# declare *new* CDS object
oligo.wt.cds <- new_cell_data_set(
  express.mat, 
  cell_metadata = cell.meta, 
  gene_metadata = gene.mat
  )

recreate.partition <- c(rep(1, length(oligo.wt.cds@colData@rownames)))
names(recreate.partition) <- oligo.wt.cds@colData@rownames
recreate.partition <- as.factor(recreate.partition)

oligo.wt.cds@clusters@listData[["UMAP"]][["partitions"]] <- recreate.partition


# (ii) cluster cells (based on seurat UMAP)
# map over additional seurat info to new CDS object
list_cluster <- seu.obj.cc.oligo.wt@active.ident
names(list_cluster) <- seu.obj.cc.oligo.wt@assays[["RNA"]]@data@Dimnames[[2]]

oligo.wt.cds@clusters@listData[["UMAP"]][["clusters"]] <- list_cluster
oligo.wt.cds@clusters@listData[["UMAP"]][["louvain_res"]] <- "NA"
oligo.wt.cds@int_colData@listData$reducedDims@listData[["UMAP"]] <- seu.obj.cc.oligo.wt@reductions[["umap"]]@cell.embeddings

# (iii) learn trajectory graph
oligo.wt.cds <- learn_graph(oligo.wt.cds, use_partition = FALSE)

pre_order.plot <- plot_cells(
  oligo.wt.cds,
  color_cells_by = 'cluster',
  cell_size = 0.6,
  label_groups_by_cluster = TRUE,
  group_label_size = 5,
  label_branch_points = FALSE,
  label_roots = FALSE
) +
  theme(legend.position = "right") +
  scale_y_continuous(limits = c(ymin, ymax))

# (iv) order cells by pseudotime
oligo.wt.cds <- order_cells(
  oligo.wt.cds,
  reduction_method = 'UMAP',
  root_cells = colnames(oligo.wt.cds[,clusters(oligo.wt.cds) == '0'])
  )

post_order.plot <- plot_cells(
  oligo.wt.cds,
  color_cells_by = 'pseudotime',
  cell_size = 0.6,
  group_label_size = 10,
  label_groups_by_cluster = FALSE,
  label_branch_points = FALSE,
  label_leaves = FALSE,
  label_roots = FALSE
) +
  theme(legend.position = "right") +
  scale_y_continuous(limits = c(ymin, ymax))

# (v) genes that change as function of pseudotime
# visualizing pseudotime in seurat
seu.obj.cc.oligo.wt$pseudotime <- pseudotime(oligo.wt.cds)
Idents(seu.obj.cc.oligo.wt) <- seu.obj.cc.oligo.wt$seurat_clusters
pseudo_feature_plot <- FeaturePlot(
  seu.obj.cc.oligo.wt, features = "pseudotime", label = T, pt.size = 1) +
  scale_y_continuous(limits = c(ymin, ymax)) +
  scale_x_continuous(limits = c(xmin, xmax))
```

### Cell Cycle Plot (WT) [Request #5]
```{r}
# Order factor levels (G1, S, G2M) to align with desired colors
phase_colors <- brewer.pal(3, "Set1")
seu.obj.cc.oligo.wt@meta.data$Phase <- factor(
  seu.obj.cc.oligo.wt@meta.data$Phase,
  levels = c("S", "G1", "G2M")
)
all_phases <- levels(seu.obj.cc.oligo.wt@meta.data$Phase)

# Check for length misalignment
if(length(all_phases) > length(phase_colors)){
  stop("Number of phases exceeds length of color palette. Choose larger palette or adjust accordingly.")
}

# Fix colors for cell-cycle phases (G1=blue, S=red, G2M=green)
names(phase_colors) <- all_phases

# Cell cycle plot
Idents(seu.obj.cc.oligo.wt) <- "Phase"
seu.obj.cc.oligo.wt.plot <- DimPlot(
  seu.obj.cc.oligo.wt, 
  label = FALSE, 
  label.size = 3,
  cols = alpha(phase_colors, 0.6),
  group.by = "Phase"
  ) + 
  ggtitle("CC-Oligodendrocytes-WT-UMAP") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-15,15))

seu.obj.cc.oligo.wt.plot
```

## Cell-Cycle Stats
```{r}
# Function cc_count() takes input seurat object and outputs matrix of the format:
#         G1  G2M   S
# WT_0d   --   --   --
# WT_4d	  --   --   --
# WT_14d  --   --   --

cc_count <- function(anon_seu_obj){

  # WT Day 0
  wt_d0_g1 <- anon_seu_obj %>%
    dplyr::filter(treatment == 'WT', time.point == '0d', Phase == 'G1') %>%
    summarise(count = n())
  
  wt_d0_g2 <- anon_seu_obj %>%
    dplyr::filter(treatment == 'WT', time.point == '0d', Phase == 'G2M') %>%
    summarise(count = n())
  
  wt_d0_s <- anon_seu_obj %>%
    dplyr::filter(treatment == 'WT', time.point == '0d', Phase == 'S') %>%
    summarise(count = n())
  

  # WT Day 4
  wt_d4_g1 <- anon_seu_obj %>%
    dplyr::filter(treatment == 'WT', time.point == '4d', Phase == 'G1') %>%
    summarise(count = n())
  
  wt_d4_g2 <- anon_seu_obj %>%
    dplyr::filter(treatment == 'WT', time.point == '4d', Phase == 'G2M') %>%
    summarise(count = n())
  
  wt_d4_s <- anon_seu_obj %>%
    dplyr::filter(treatment == 'WT', time.point == '4d', Phase == 'S') %>%
    summarise(count = n())
  
  # WT Day 14
  wt_d14_g1 <- anon_seu_obj %>%
    dplyr::filter(treatment == 'WT', time.point == '14d', Phase == 'G1') %>%
    summarise(count = n())
  
  wt_d14_g2 <- anon_seu_obj %>%
    dplyr::filter(treatment == 'WT', time.point == '14d', Phase == 'G2M') %>%
    summarise(count = n())
  
  wt_d14_s <- anon_seu_obj %>%
    dplyr::filter(treatment == 'WT', time.point == '14d', Phase == 'S') %>%
    summarise(count = n())

  cc_table <- data.frame(
    G1 = c(wt_d0_g1$count[1], 
           wt_d4_g1$count[1],
           wt_d14_g1$count[1]),
    G2M = c(wt_d0_g2$count[1], 
           wt_d4_g2$count[1],
           wt_d14_g2$count[1]),
    S = c(wt_d0_s$count[1], 
           wt_d4_s$count[1],
           wt_d14_s$count[1]))
  
  return(cc_table)
}

# Function cc_percent() takes input matrix of format equivalent to cc_count() and outputs matrix of same format but with respective percentage values in cell
cc_percent <- function(cc_count_table){
  
  wt_0d_count <- cc_count_table[1,]
  wt_4d_count <- cc_count_table[2,]
  wt_14d_count <- cc_count_table[3,]
  
  wt_0d_perc <- (wt_0d_count / sum(wt_0d_count)) * 100
  wt_4d_perc <- (wt_4d_count / sum(wt_4d_count)) * 100
  wt_14d_perc <- (wt_14d_count / sum(wt_14d_count)) * 100


  cc_perc_table <- rbind(wt_0d_perc, 
                         wt_4d_perc, 
                         wt_14d_perc)
  return(cc_perc_table)
}
```

```{r}
row_names <- c("WT_0d", "WT_4d", "WT_14d")

# Microglia
# Create count matrix table
mg_cc_table <- cc_count(seu.obj.cc.mg@meta.data)
row.names(mg_cc_table) <- row_names
mg_cc_table

# Confirm total count observations in cc_table 
# matches total observations in seurat object
match_status <- sum(unlist(mg_cc_table)) == dim(seu.obj.cc.mg@meta.data)[1]
cat("Does cc_table total obs equal total obs of seurat object?:", match_status)

# Create percent matrix table
mg_cc_perc_table <- cc_percent(mg_cc_table)
```


```{r}
# Astrocytes
row_names <- c("WT_0d", "WT_4d", "WT_14d")

# Create count matrix table
astro_cc_table <- cc_count(seu.obj.cc.astro@meta.data)
row.names(astro_cc_table) <- row_names
astro_cc_table


# Confirm total count observations in cc_table 
# matches total observations in seurat object
match_status <- sum(unlist(astro_cc_table)) == dim(seu.obj.cc.astro@meta.data)[1]
cat("Does cc_table total obs equal total obs of seurat object?:", match_status)

# Create percent matrix table
astro_cc_perc_table <- cc_percent(astro_cc_table)
```


```{r}
# Oligodendrocytes
row_names <- c("WT_0d", "WT_4d", "WT_14d")

# Create count matrix table
oligo_cc_table <- cc_count(seu.obj.cc.oligo@meta.data)
row.names(oligo_cc_table) <- row_names
oligo_cc_table


# Confirm total count observations in cc_table 
# matches total observations in seurat object
match_status <- sum(unlist(oligo_cc_table)) == dim(seu.obj.cc.oligo@meta.data)[1]
cat("Does cc_table total obs equal total obs of seurat object?:", match_status)

# Create percent matrix table
oligo_cc_perc_table <- cc_percent(oligo_cc_table)
```

# Selected gene expression plot for Microglia WT
```{r, eval = FALSE}
# subset out Microglia
Idents(seu.cc.sub.wt) <- "celltype"

wt.mg <- subset(seu.cc.sub.wt, idents="Microglia")

# wt
Spp1.wt.feat <- FeaturePlot(wt.mg, features = "Spp1") + 
  theme(legend.position = "none")

Axl.wt.feat <- FeaturePlot(wt.mg, features = "Axl") + 
  theme(legend.position = "none")

Itgax.wt.feat <- FeaturePlot(wt.mg, features = "Itgax") + 
  theme(legend.position = "none")

Apoe.wt.feat <- FeaturePlot(wt.mg, features = "Apoe") + 
  theme(legend.position = "none")

Tmem119.wt.feat <- FeaturePlot(wt.mg, features = "Tmem119") + 
  theme(legend.position = "none")

Trem2.wt.feat <- FeaturePlot(wt.mg, features = "Trem2") + 
  theme(legend.position = "none")


# plot grid total
feat.grid.mg <- grid.arrange(
  textGrob("WT"), Spp1.wt.feat, Axl.wt.feat, Itgax.wt.feat, Apoe.wt.feat, Tmem119.wt.feat, Trem2.wt.feat,
  ncol=7, top = "Microglia Select Genes")

```


# Selected gene expression plot for Astrocytes WT
```{r}
# subset out Astrocytes
seu.obj.cc.wt.astro <- subset(seu.obj.cc.wt, idents="Astrocytes")

# wt
Aqp4.wt.feat <- FeaturePlot(seu.obj.cc.wt.astro, features = "Aqp4") + 
  theme(legend.position = "none")

Gfap.wt.feat <- FeaturePlot(seu.obj.cc.wt.astro, features = "Gfap") + 
  theme(legend.position = "none")

Gja1.wt.feat <- FeaturePlot(seu.obj.cc.wt.astro, features = "Gja1") + 
  theme(legend.position = "none")

Slc1a2.wt.feat <- FeaturePlot(seu.obj.cc.wt.astro, features = "Slc1a2") + 
  theme(legend.position = "none")

Aldh1l1.wt.feat <- FeaturePlot(seu.obj.cc.wt.astro, features = "Aldh1l1") + 
  theme(legend.position = "none")



# plot grid total
grid.arrange(textGrob("WT"), Aqp4.wt.feat, Gfap.wt.feat, Gja1.wt.feat, Slc1a2.wt.feat, Aldh1l1.wt.feat,
             ncol=6, top = "Astrocytes Select Genes")

```

# Panel Plotting - via Gene
Create separate dataframe (meta.df) to manipulate and use with ggplot
```{r}
# metadata as new dataframe
meta.df <- seu.obj.cc.astro.wt@meta.data

# Make cell_IDs a column in Seurat Obj meta data table
meta.df$cell_ID <- rownames(seu.obj.cc.astro.wt@meta.data)
```

```{r}
meta.df <- 
  seu.obj.cc.astro.wt@reductions$umap@cell.embeddings %>%
  as.data.frame() %>%
  tibble::rownames_to_column(var='barcode') %>%
  inner_join(meta.df, by=c('barcode'='cell_ID')) %>%
  dplyr::rename("cell_ID" = "barcode")
```

```{r}
# convert count matrix from type dgcMatrix to type data.frame
rna.count.mat <- as.data.frame(seu.obj.cc.astro.wt@assays$RNA@counts)


# Transpose matrix so that Gene Expressions become Column Names, 
# Cell ID's become rownames, and Cell ID's become a column for identification
roti.rna.count.mat <- as.data.frame(t(rna.count.mat))
roti.rna.count.mat$cell_ID <- rownames(roti.rna.count.mat)
```

### Append Gene Count Function
```{r}
append_gene_count <- function(met.df, main.matrix, gene.list){
  
  gene.vec <- vector("numeric", length=0)
  
  for(gene in gene.list){
    gene.vec <- main.matrix %>%
      dplyr::select(cell_ID, !!gene)
  
    
    met.df <- dplyr::right_join(met.df, gene.vec, by = "cell_ID")
  }
  
  return(met.df)
}
```

Function checks whether specific gene can be found
```{r}
check_gene <- function(main.matrix, gene.list){
  for(gene in gene.list){
    search.result <- gene %in% colnames(main.matrix)
    if(search.result == TRUE){
      cat(gene, "\t","\t", "status in matrix is:", search.result, "\n")
    } else {
      cat(gene, "\t","\t","\t", "status in matrix is:", search.result, "\n")
    }
  }
}


```
## Append target genes to metadata
T Cells         Cd3e
B Cells         Cd19
Neutrophils     Ly6g, S100a9
Mono/Mac        Adgre1, Cd86, Cd40  
Microglia       Tmem119, P2ry12, Gpr34, Sall1  
Astrocytes      Gfap, Gja1, Aqp4, Slc1a2, Aldh1l1, Rax
Schwann Cells   Npy, Gja1, Cldn5, Rax
Pericytes       Kcnj8, Cspg4
Endothelial Cells     Cldn5, Rax
VSMC            Acta2, Cspg4
VLMC            Slc6a13, Pdgfra, Gja1
NRP             Sox11, Cdk1
ImmN            Sox11, Syt1
TNC             Rax, Gja1
Oligodendro.    Pdgfra, Cldn11, Cspg4, Syt1, Rax
CPC             Folr1, Ttr, Pdgfra, Syt1, Rax
```{r}
target.genes <- c("Cd3e",
                  "Cd19",
                  "Ly6g", "S100a9",
                  "Adgre1", "Cd86", "Cd40",
                  "Tmem119", "P2ry12", "Gpr34", "Sall1",
                  "Gfap", "Gja1", "Aqp4", "Slc1a2", "Aldh1l1", "Rax",
                  "Npy", "Cldn5",
                  "Kcnj8", "Cspg4",
                  "Cldn5",
                  "Acta2",
                  "Slc6a13", "Pdgfra",
                  "Sox11", "Cdk1",
                  "Syt1",
                  "Cldn11",
                  "Folr1", "Ttr",
                  "Rbfox3", "Tubb3", "Map2"
                  )
#check_gene(roti.rna.count.mat, target.genes)

meta.df <- append_gene_count(meta.df, 
                                 roti.rna.count.mat, 
                                 target.genes
                                 )

# genes that were not detected in the matrix are listed in the bench vector below:
# benched.genes <- c("Cspd4", "NeuN", "Tuj1")
```

Update column variables to appropriate datatypes for ggplot graphing
```{r}
meta.df <- 
  meta.df %>%
  mutate(across(Cd3e:Map2,
                ~ as.numeric(.x)
                )
  ) %>%
  mutate(time.point = factor(time.point,
                            levels=c("0d", "4d", "14d")
                            )
        )
```


meta.df plotted via ggplot
Plotting w/ggplot: Microglia genes (Tmem119, P2ry12, Gpr34, Sall1)
Astrocyte genes: ("Gfap", "Gja1", "Aqp4", "Slc1a2", "Aldh1l1")
```{r}
meta.df %>%
  ggplot(aes(x = UMAP_1,
            y = UMAP_2,
            alpha = 0.2,
            color = log(Gfap)
          )
  ) + 
  geom_point(size = 1.2,
            stroke = 0, 
            shape = 16) +
  scale_color_gradient(low="#f5f1c4", high="#f20a0a") +
  theme_bw() +
  coord_fixed(ratio = 1,  xlim = c(-10,20), ylim = c(-20,10)) + 
  facet_grid(treatment ~ time.point)
```

```{r}
meta.df %>%
  ggplot(aes(x = UMAP_1,
            y = UMAP_2,
            alpha = 0.2,
            color = log(Slc1a2)
          )
  ) + 
  geom_point(size = 1.2,
            stroke = 0, 
            shape = 16) +
  scale_color_gradient(low="#f5f1c4", high="#f20a0a") +
  theme_bw() +
  coord_fixed(ratio = 1, xlim = c(-20,20), ylim = c(-20,20)) + 
  facet_grid(treatment ~ time.point)
```

# [Request #4] Panel Plotting - via TimePoint

## Microglia
```{r}
# subset mg seu obuject to WT
Idents(seu.obj.cc.mg) <- "treatment"
seu.obj.cc.mg.wt <- subset(seu.obj.cc.mg, idents="WT")

# lock colors for seurat clusters for consistent panel plotting
all_clusters <- levels(seu.obj.cc.mg.wt@meta.data$seurat_clusters)
cluster_colors <- c(
  "Population_A" = "#1976d2",
  "Population_B" = "#66bb6a",
  "Population_C" = "#9C27B0",
  "Population_D" = "orange")
#cluster_colors <- brewer.pal(4, "RdYlBu")
if(length(all_clusters) > length(cluster_colors)){
  stop("Number of clusters exceeds the length of the color palette. Choose a larger palette or adjust accordingly.")
}
names(cluster_colors) <- all_clusters
```

### WT
```{r}
# merged (0d, 4d, 14d) plot
mg.merge <- DimPlot(
  seu.obj.cc.mg.wt,
  pt.size = .8,
  cols = alpha(cluster_colors, 0.3),
  group.by = "seurat_clusters"
  ) + 
  ggtitle("Merged") +
  # theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-10,10))

# extract legend as own grob to be attached at end panel
mg.merge.grob <- ggplotGrob(mg.merge)
legend.grob <- gtable::gtable_filter(mg.merge.grob, "guide-box")
mg.merge <- mg.merge + theme(legend.position = "none")


Idents(seu.obj.cc.mg.wt) <- "time.point"
mg.wt.0d <- subset(seu.obj.cc.mg.wt, idents="0d")
mg.wt.4d <- subset(seu.obj.cc.mg.wt, idents="4d")
mg.wt.14d <- subset(seu.obj.cc.mg.wt, idents="14d")

# 0d plot
mg.wt.0d.plot <- DimPlot(
  mg.wt.0d,
  pt.size = 0.8,
  cols = alpha(cluster_colors, 0.5),
  group.by = "seurat_clusters"
  ) + 
  ggtitle("0 days") +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-10,10))

# 4d plot
mg.wt.4d.plot <- DimPlot(
  mg.wt.4d,
  pt.size = 0.8,
  cols = alpha(cluster_colors, 0.3),
  group.by = "seurat_clusters"
  ) + 
  ggtitle("4 days") +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-10,10))

# 14d plot
mg.wt.14d.plot <- DimPlot(
  mg.wt.14d,
  pt.size = 0.8,
  cols = alpha(cluster_colors, 0.3), 
  group.by = "seurat_clusters"
  ) + 
  ggtitle("14 days") +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-10,10))

mg.label <- grid::textGrob(
  "Microglia (WT)",
  gp=grid::gpar(fontsize=20, fontface="bold"),
  rot = 90
  )

grid.arrange(
    arrangeGrob(
        mg.merge, 
        mg.wt.0d.plot, 
        mg.wt.4d.plot, 
        mg.wt.14d.plot,
        ncol = 4,
        left = mg.label,
        right = legend.grob))


```

## Astrocytes
```{r}
# subset astro seu object to WT
Idents(seu.obj.cc.astro) <- "treatment"
seu.obj.cc.astro.wt <- subset(seu.obj.cc.astro, idents="WT")

# lock colors for seurat clusters for consistent panel plotting
all_clusters <- levels(seu.obj.cc.astro.wt@meta.data$seurat_clusters)
cluster_colors <- brewer.pal(7, "Dark2")
if(length(all_clusters) > length(cluster_colors)) {
  stop("Number of clusters exceeds the length of the color palette. Choose a larger palette or adjust accordingly.")
}
names(cluster_colors) <- all_clusters
```

### WT
```{r}
# merged  (0d, 4d, 14d) plot
astro.merge <- DimPlot(
  seu.obj.cc.astro.wt,
  pt.size = 0.3,
  cols = alpha(cluster_colors, 0.3),
  group.by = "seurat_clusters"
  ) + 
  ggtitle("Merged") +
  #theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-20,10))

# extract legend as own grob to be attached at end panel
astro.merge.grob <- ggplotGrob(astro.merge)
legend.grob <- gtable::gtable_filter(astro.merge.grob, "guide-box")
astro.merge <- astro.merge + theme(legend.position = "none")

Idents(seu.obj.cc.astro.wt) <- "time.point"
astro.wt.0d <- subset(seu.obj.cc.astro.wt, idents="0d")
astro.wt.4d <- subset(seu.obj.cc.astro.wt, idents="4d")
astro.wt.14d <- subset(seu.obj.cc.astro.wt, idents="14d")


# 0d plot
astro.wt.0d.plot <- DimPlot(
  astro.wt.0d, 
  pt.size = 0.3,
  cols = alpha(cluster_colors, 0.3),
  group.by = "seurat_clusters"
  ) + 
  ggtitle("0 days") +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-20,10))

# 4d plot
astro.wt.4d.plot <- DimPlot(
  astro.wt.4d,
  pt.size = 0.3,
  cols = alpha(cluster_colors, 0.3),
  group.by = "seurat_clusters"
  ) + 
  ggtitle("4 days") +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-20,10))

# 14d plot
astro.wt.14d.plot <- DimPlot(
  astro.wt.14d, 
  pt.size = 0.3,
  cols = alpha(cluster_colors, 0.3),
  group.by = "seurat_clusters"
  ) + 
  ggtitle("14 days") +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-20,10))

astro.label <- grid::textGrob(
  "Astrocytes (WT)",
  gp=grid::gpar(fontsize=12, fontface="bold"),
  rot = 90
  )


grid.arrange(
    arrangeGrob(
        astro.merge, 
        astro.wt.0d.plot, 
        astro.wt.4d.plot, 
        astro.wt.14d.plot,
        ncol = 4,
        left = astro.label,
        right = legend.grob))

```

## Oligodendrocytes
```{r}
Idents(seu.obj.cc.oligo) <- "treatment"
seu.obj.cc.oligo.wt <- subset(seu.obj.cc.oligo, idents="WT")

# lock colors for seurat clusters for consistent panel plotting
all_clusters <- levels(seu.obj.cc.oligo.wt@meta.data$seurat_clusters)
cluster_colors <- brewer.pal(5, "Dark2")
if(length(all_clusters) > length(cluster_colors)) {
  stop("Number of clusters exceeds the length of the color palette. Choose a larger palette or adjust accordingly.")
}
names(cluster_colors) <- all_clusters
```

### WT
```{r}
# oligo merge plot
oligo.merge <- DimPlot(
  seu.obj.cc.oligo.wt, 
  pt.size = 0.4,
  group.by = "seurat_clusters",
  cols = alpha(cluster_colors, 0.4)) + 
  ggtitle("Merged") +
  #theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-15,15))

# extract legend as own grob to be attached at end panel
oligo.merge.grob <- ggplotGrob(oligo.merge)
legend.grob <- gtable::gtable_filter(oligo.merge.grob, "guide-box")
oligo.merge <- oligo.merge + theme(legend.position = "none")

Idents(seu.obj.cc.oligo.wt) <- "time.point"
oligo.wt.0d <- subset(seu.obj.cc.oligo.wt, idents="0d")
oligo.wt.4d <- subset(seu.obj.cc.oligo.wt, idents="4d")
oligo.wt.14d <- subset(seu.obj.cc.oligo.wt, idents="14d")

# 0d plot
oligo.wt.0d.plot <- DimPlot(
  oligo.wt.0d, pt.size = 0.4,
  cols = alpha(cluster_colors, 0.4),
  group.by = "seurat_clusters"
  ) + 
  ggtitle("0 days") +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-15,15))
  
# 4d plot
oligo.wt.4d.plot <- DimPlot(
  oligo.wt.4d,pt.size = 0.4,
  cols = alpha(cluster_colors, 0.4),
  group.by = "seurat_clusters"
  ) + 
  ggtitle("4 days") +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-15,15))

# 14d plot
oligo.wt.14d.plot <- DimPlot(
  oligo.wt.14d, pt.size = 0.4,
  cols = alpha(cluster_colors, 0.4),
  group.by = "seurat_clusters"
  ) + 
  ggtitle("14 days") +
  theme(legend.position = "none") +
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-10,20), ylim = c(-15,15))

oligo.label <- grid::textGrob(
  "Oligodendrocytes (WT)",
  gp=grid::gpar(fontsize=12, fontface="bold"),
  rot = 90
  )


grid.arrange(
    arrangeGrob(
        oligo.merge, 
        oligo.wt.0d.plot, 
        oligo.wt.4d.plot, 
        oligo.wt.14d.plot,
        ncol = 4,
        left = oligo.label,
        right = legend.grob
    )
)

```

Issue: want to plot 0d, 4d, 14d side by side with 
Custom plotting with seu object metadata.

Extra code for original panel plotting (no y-axis)
```{r}
# enables side to side plotting (not used since Y-axis wanted)
time.list <- list("14d", "4d", "0d")
mg.0thru14 <- DimPlot(seu.obj.cc.mg.wt, 
                 split.by = 'ident',
                 order = time.list
                 ) + 
  theme(plot.title = element_text(hjust = 0.5)) +
  coord_fixed(ratio = 1, xlim = c(-15,10), ylim = c(-15,10))

grid.arrange(mg.merge, mg.0thru14, ncol = 2)
```

# [Request #7] Differential Gene Expression
*0d WT microglia vs 4d WT microglia*
*0d WT microglia vs 14d WT microglia*
```{r}
Idents(seu.cc.sub.wt) <- "celltype"
seu.cc.sub.wt.mg <- subset(seu.cc.sub.wt, idents = 'Microglia')

Idents(seu.cc.sub.wt.mg) <- "time.point"
wt.mg.markers.4d.0d <- FindMarkers(seu.cc.sub.wt.mg, ident.1 = "4d", ident.2 = "0d")

wt.mg.markers.14d.0d <- FindMarkers(seu.cc.sub.wt.mg, ident.1 = "14d", ident.2 = "0d")
```

*0d WT astrocyte vs 4d WT astrocyte*
*0d WT astrocyte vs 14d WT astrocyte*
```{r}
Idents(seu.cc.sub.wt) <- "celltype"
seu.cc.sub.wt.astro <- subset(seu.cc.sub.wt, idents = 'Astrocytes')

Idents(seu.cc.sub.wt.astro) <- "time.point"
wt.astro.markers.4d.0d <- FindMarkers(seu.cc.sub.wt.astro, ident.1 = "4d", ident.2 = "0d")

wt.astro.markers.14d.0d <- FindMarkers(seu.cc.sub.wt.astro, ident.1 = "14d", ident.2 = "0d")
```

*0d WT oligodendrocyte vs 4d WT oligodendrocyte*
*0d WT oligodendrocyte vs 14d WT oligodendrocyte*
```{r}
Idents(seu.cc.sub.wt) <- "celltype"
seu.cc.sub.wt.oligo <- subset(seu.cc.sub.wt, idents = 'Oligodendrocytes')

Idents(seu.cc.sub.wt.oligo) <- "time.point"
wt.oligo.markers.4d.0d <- FindMarkers(seu.cc.sub.wt.oligo, ident.1 = "4d", ident.2 = "0d")

wt.oligo.markers.14d.0d <- FindMarkers(seu.cc.sub.wt.oligo, ident.1 = "14d", ident.2 = "0d")
```


```{r}
dir_path <- "/hpc/home/jcy13/scRNA_crypto_Jae/Datasets"

# Put the data tables into a named list
table.list <- list(mg.markers.4d.0d = wt.mg.markers.4d.0d,
                   mg.markers.14d.0d = wt.mg.markers.14d.0d,
                   astro.markers.4d.0d = wt.astro.markers.4d.0d,
                   astro.markers.14d.0d = wt.astro.markers.14d.0d,
                   oligo.markers.4d.0d = wt.oligo.markers.4d.0d,
                   oligo.markers.14d.0d = wt.oligo.markers.14d.0d)

# Loop through the list and save each data table as a CSV
for (name in names(table.list)) {
  # Construct the file name based on the list name
  file_path <- file.path(dir_path, paste0(name, ".csv"))
  
  # Save the data table as a CSV
  fwrite(table.list[[name]], file = file_path, row.names = TRUE)
}
```
